Layered Tin Chalcogenides SnS and SnSe: Lattice Thermal Conductivity Benchmarks and Thermoelectric Figure of Merit

Tin sulfide (SnS) and tin selenide (SnSe) are attractive materials for thermoelectric conversion applications. Favorable small band gap, high carrier mobility, large Seebeck coefficient, and remarkably low lattice thermal conductivity are a consequence of their anisotropic crystal structure of symmetry Pnma, made of corrugated, black phosphorus-like layers. Their internal lattice dynamics combined with chemical bond softening in going from SnS to SnSe make for subtle effects on lattice thermal conductivity. Reliable prediction of phonon transport for these materials must therefore include many-body effects. Using first principles methods and a transferable tight-binding potential for frozen phonon calculations, here, we investigate the evolution of thermal lattice conductivity and thermoelectric figure of merit in Pnma-SnS and -SnSe, also including the high-temperature Cmcm-SnS phase. We show how thermal conductivity lowering in SnS at higher temperatures is largely due to dynamic phonon softening ahead of the Pnma–Cmcm structural phase transition. SnS becomes more similar to SnSe in its lifetime and mean free path profiles as it approaches its high-temperature Cmcm phase. The latter nonetheless intrinsically constraints phonon group velocity modules, preventing SnS to overtake SnSe. Our analysis provides important insights and computational benchmarks for optimization of thermoelectric materials via a more efficient computational strategy compared to previous ab initio attempts, one that can be easily transferred to larger systems for further thermoelectric materials nanoengineering. The good description of anharmonicity at higher temperatures inherent to the tight-binding potential yields calculated lattice conductivity values that are in very good agreement with experiments.


■ INTRODUCTION
The increased energy demand and migration away from fossil fuel combustion call for innovative material solutions for sustainable development. Automotive exhaust and industrial processes waste about two-thirds of the total energy into heat. Furthermore, harvesting heat from battery packs, fuel, and photovoltaic cells increases the overall device energy conversion efficiency and process sustainability and circularity. Finding high-performance thermoelectric materials capable of directly and reversibly converting heat to electrical energy is therefore a task of top priority. 1−7 The efficient control of thermoelectric energy conversion processes entails the ability to assemble materials with tailored thermal transport properties. 8 The thermoelectric efficiency at a given temperature is expressed by the dimensionless figure of merit ZT = S 2 σT/κ. Here, S is the Seebeck coefficient, σ is the electrical conductivity, and κ is the thermal conductivity, which is the sum of lattice and electronic contributions κ = κ latt + κ e . This expression suggests a large Seebeck coefficient and high electrical conductivity as requirements for good thermoelectric materials, as well as low thermal conductivity. Nonetheless, these conditions cannot be immediately translated into a material design strategy due to the interdependence of the physical parameters contributing to ZT.
With lead chalcogenides firmly established in the literature as strong thermoelectrics, 6,9,10 tin chalcogenides have emerged as suitable replacements exhibiting lower toxicity than Pbbased materials. Binary tin chalcogenides SnX (X = S, Se) owe their versatility and property tunability to their distorted NaCl geometry, which yields an anisotropic layered structure of the Pnma symmetry. With a small band gap, high carrier mobility, a large Seebeck coefficient, and extremely low lattice thermal conductivity, SnSe has been the focus of an intense search for high-figure-of-merit materials. 11 With values as low as 0.25 W m −1 K −1 above 800 K, SnSe outperforms top composite thermoelectric materials, without any need for grain engineering, 12 nanoinclusions, or nanostructuring. Morphology, vacancies, and synthetic approaches can sensibly affect thermal conductivity. 13,14 At high temperatures, both Pnma-SnS and -SnSe undergo a structural phase transition to a Cmcm structure, in which the c-axis shortens to more closely resemble the b-axis. This follows from the compression of the intrinsically deformable, zig-zag "accordion-like" layers' geometry and from antiparallel layer displacements. Interestingly, it is for temperatures near to and above the phase transition threshold that the highest ZT figures are measured. While low lattice thermal conductivities can be expected from the anisotropic, layered structure alone, 15 the relevance of several structures toward enhancing thermoelectricity is therefore very intriguing. On crossing the phase boundary, favorable electronic properties are maintained, while layer compression and reshuffling further assist phonon softening. 15 Detailed diffraction studies suggest simultaneous positive and negative Seebeck effects along different directions as a consequence of structural changes induced by pressure. 16 Low thermal transport in Cmcm-SnS is associated with phonon softening and enhanced three-phonon scattering, affecting in particular low-frequency modes. 17 The lighter homologue SnS has a comparatively smaller thermoelectric figure of merit at lower temperatures, while ab initio calculations suggest the Cmcm phase to be competitive in SnSe at higher temperatures, 18 confirming the subtle structure−property relationship at work in this class of materials. 19 Pressure strongly affects transport in SnS, inducing superconductivity, an effect already noticed in black phosphorus and associated with an electronic fingerprint of a good thermoelectric response. 20,21 Heat transport properties of SnS and SnSe, individually or comparatively, have been calculated from first principles in several contexts. 15,17−19,22−27 The interplay of the anisotropic structure, temperature-induced phonon softening, and manybody phonon scattering has been indicated as an important factor in explaining remarkably low lattice conductivities. 15,17,25 While both Pnma-SnS and -SnSe are consistently evaluated as poor thermal conductors, their relative ranking has been more challenging to determine from first principles calculations. 23 In the search for top-performing thermoelectrics, complex materials are indicated as technology relevant, comprising nanoincluded composites, mesoscale grain boundaries, or superlattice engineering to effectively act as phonon scatterers. 8,10 To further improve on the figure of merit of tin chalcogenides, it is therefore desirable to have a reliable yet computationally efficient baseline for thermal transport benchmark calculations, one that can be easily implemented on a larger structural scale, for which ab initio approaches are impractical. As with all thermoelectric materials, the question of doping is also pertinent, and obtaining sufficient doping levels without detriment to the lattice dynamics is a delicate task. Dopants commonly used are sodium 28 and cadmium 13 for p-type and chlorine 29 for n-type.
In this work, we set such a baseline by electronic and thermal transport in Pnma-SnS and -SnSe, including a comparison with Cmcm-SnS. The electronic part of the thermoelectric figure of merit is obtained from the Boltzmann transport equation based on the Wannier function (WF) interpolation of the electronic bands. 30 Phonon frequencies are calculated from a transferable tight-binding interatomic potential, 31 while lattice thermal conductivity is obtained from the linearized phonon Boltzmann equation (LBTE) using a supercell approach. 32,33 Thermal transport anisotropy and ranking are obtained in good agreement with experiments, while the increasing similarity of SnS and SnSe is mapped onto a common Cmcm structural motif, which becomes dynamically important at high temperatures. ■ METHODS Electronic Structure and Phonons. Bloch states were calculated using density functional theory in the GGA-PBE approximation, using Vanderbilt ultrasoft Sn, S, and Se standard solid-state pseudopotentials (SSSP, 34 efficiency branch) available from the materials cloud (materialscloud.org). Pseudopotentials included 5s 2 5p 2 4d 10 for Sn and 3s 2 3p 2 for S and Se. The plane-wave expansion of valence electron wave functions and charge density used kinetic energy cutoffs of 70 and 560 Ry, respectively. Calculations were performed using the pw.x code of the Quantum Espresso (QE) package v. 6.7. 35 SCF calculations were run on a k-grid of 3 × 8 × 8 data points. Cell parameters and atomic coordinates were relaxed below 1 × 10 −8 Ry and 1 × 10 −6 Ry/Bohr in energy and forces, respectively. Bloch states and overlaps for maximally localized WF (MLWF) calculations were evaluated on a Monkhorst− Pack 4 × 11 × 10 k-point mesh for both SnS and SnSe.
A tight-binding band structure model of Bloch functions was obtained from interpolating MLWFs, calculated with the Wannier90 package. 36 As an initial guess for Wannier functions, p orbitals were used for Sn, S, and Se (24 MLWFs), and s projectors were included for Sn, excluding energetically low lying 3s bands (28 MLWF model). Iterative spread minimization provided real-valued, maximally localized WFs using a convergence tolerance of 1.0 × 10 −10 Ω.
The Seebeck tensor S depends on the doping level of the material through the chemical potential μ and the temperature T. Collisions in the semiclassical Boltzmann transport equation are mapped into band n and k vector-specific relaxation times, τ nk . In the single relaxation approximation, a constant relaxation time is assumed for every band and momentum, τ = τ nk .
The postw90 package of Wannier90 was used to interpolate band velocities v and calculate the transport distribution f u n c t i o n Electrical conductivity σ, Seebeck coefficient S, and electronic thermal conductivity κ = K − σS 2 T tensors were calculated semiclassically with the Boltzmann transport equation, using the BoltzWann module of Wannier90. 30 The dependency on temperature and chemical potential is given as Thermal Conductivity. The linearized phonon Boltzmann transport equation (LBTE) was solved using the single-mode relaxation approximation (SMRT/RTA). 32,33 Therein, the lattice thermal conductivity contributed by phonon modes was calculated from mode heat capacity, phonon group v e l o c i t y , a n d r e l a x a t i o n t i m e a c c o r d i The mode heat capacity C λ for phonon frequency ω λ is given as follows The many-body nature of the scattering process introduces three-phonon interaction contributions to phonon mode lifetimes. The corresponding three-phonon interaction strengths can be computed from phonon harmonic frequencies, displacement eigenvectors, and third-order force constant matrices. 32,33 The imaginary part of the self-energy, computed from three-body phonon scattering, takes the form n n n n Φ λλ′λ″ is the interaction strength between three scattering λ, λ′, and λ″ phonons, while n λ is the temperature-dependent equilibrium phonon occupation number, which is given as follows: relates to the phonon lifetime as , where 2 ( ) is the phonon linewidth of the λ phonon mode. For the phonon and lattice thermal conductivity calculations, PHONOPY 39 and PHONO3PY 32 software packages were employed.
In the SMRT/RTA approach, it is assumed that phonon relaxation can be replaced by phonon lifetime. This approximation can be checked by performing a computational demanding but more accurate full solution of the LBTE, according to ref 40. This was done for both SnS and SnSe (Pnma) to assess the numerical quality of the SMRT/RTA model for this class of compounds (see Tables S1 and S2 in the Supporting Information, SI).
Second-and third-order force constants were calculated based on the supercell approach with finite atomic displacements of 0.03 Å. 32 The use of larger supercells is in general important to compute phonon−phonon scattering channels with better accuracy, as spurious imaginary modes may appear at low frequencies around Γ. For the Pnma structures, we used a 3 × 3 × 3 supercell for third-order force calculations and a 4 × 8 × 8 supercell for second-order force constants. Cmcm-SnS used a 3 × 3 × 3 supercell for third-order force calculations and a 6 × 6 × 6 supercell for the second-order force constants. Quantum Espresso phonon calculations (PHONOPY) were based on a 2 × 4 × 4 supercell, while CP2K-xTB phonons used a 3 × 3 × 3 supercell. A large supercell number is required for the calculation of third-order force constants. While this task   42 using a modified interface to PHONO3PY. Geometry relaxed crystal structures were used to construct a supercell for the thermal conductivity calculation (see Figure  S1 in the SI). Cell parameters and atomic positions were relaxed using the BFGS scheme on a k-grid of 6 × 8 × 8, including long-range contributions to the forces via Ewald summation. The electronic band structures of SnS and SnSe are shown in Figure 2. The highest valence bands below the Fermi level show a characteristic "pudding mold"-like structure (Γ−Z) in SnS, which is also present in SnSe but with a more pronounced energy offset between the two maxima. A second valence band maximum is visible in correspondence with the conduction band minimum along Y−Γ. Their cusps connect two regions of linear E(k) dependence on k, like in an "opened Dirac cone." The Z−U segment in SnS lies above the Y−Γ valence maximum cusp, which is not the case for SnSe, a feature also noticed in quasiparticle GW calculations. 25 The second lowest CB minimum for both SnS and SnSe is found at Γ. Intralayer bands are more dispersive than interlayer ones. SnS shows a direct band gap of 1.08 eV and an indirect band gap of 0.92 eV, in agreement with literature values. 43 Phonons. Harmonic phonons enter the calculation of three-phonon interaction strengths. Previous ab initio phonon calculations 25 indicate similar phonon band dispersions, with optical frequencies higher for SnS due to lighter S mass. The exceptionally low thermal conductivity of SnSe in particular is connected to optical modes softening. 11 If on the one hand the crystal axis anisotropy of thermal transport is very well accounted for by ab initio approaches, comparison with experiments has been so far less conclusive regarding absolute values and relative ranking. 23 To calculate phonons and thermal conductivity, we therefore based force calculations on a transferable, semiempirical third-order tight-binding potential, GFN-xTB 41 as implemented in CP2K. 42 Figure  S1.
Tight-binding forces bias harmonic and cubic force constants to larger lattice thermal conductivity values while maintaining the same directional ranking within each structure and relative alignment as in experiments, κ latt (SnS) > κ latt (SnSe) (see Figure S2). Pnma-SnS Quantum Espresso (QE) phonon dispersions compare in general well with the GFN-xTB results, the latter showing increased dispersion along U− R−Z (in agreement with ref 25) but narrower dispersion above the phonon band gap. However, CP2K phonons are in general The Journal of Physical Chemistry C pubs.acs.org/JPCC Article harder compared to QE and published phonon spectra. 23,25,47 Forces were therefore mollified aiming at a better match of the acoustic CP2K band dispersions with QE results (based on SnS, see Figure S3), as acoustic frequencies significantly contribute to the lattice thermal conductivity. 48,49 All xTB forces were accordingly scaled to 25% of their pristine value for all three structures considered in this work. The resulting SnS and SnSe phonon spectra are shown in     (Figure 4a). For Pnma-SnSe, the axis-resolved lattice thermal conductivities κ latt at 300 K are 0.3374, 0.7726, and 0.5652 W m −1 K −1 for the a-, b-, and c-axis, respectively, with an averaged/isotropic lattice thermal conductivity κ iso of 0.5584 W m −1 K −1 . While the axial component sequence κ b > κ c > κ a is the same order as for SnS, the anisotropy of lattice thermal transport is larger than that in SnS due to a more pronounced relative lowering of κ c ( Figure  4b). The calculated values at 300 K are in agreement with experimental values (0.46, 0.70, and 0.68 W m −1 K −1 along the a-, b-, and c-axis, 15 respectively), including their decrease to lowest values at high temperatures (<0.2 W m −1 K −1 ).
SnS displays longer optical phonon relaxation times (>1.3 THz) than SnSe, for which sizable (>10 ps) relaxation times are found at lower frequencies only, see Figure 5. Optical frequency relaxation times are in general shorter in SnSe than those in SnS. Taken together, the loss of hotspot lifetime density from SnS to SnSe (Figure 5a,b) and the compactification of lifetimes in SnSe explain differences in calculated κ latt , indicating a broader impact over frequencies in the latter and not only a localized effect around lower/lowest frequencies. Furthermore, longer lifetimes at higher frequencies explain the more pronounced temperature effect on κ a in SnS compared to SnSe, for which additional thermal conduction lowering is controlled by rapidly relaxing optical modes.
Full Figure of Merit. In the single-lifetime approximation, the chemical potential μ was chosen to match the experimental Seebeck coefficients at a given temperature followed by τ fit on experimental electrical conductivity. In combination with lattice thermal conductivity, the full thermoelectric figure of merit ZT in SnS and SnSe can be evaluated and compared at different temperatures.
The relaxation times at 300 K τ = 24.9394 fs (SnS) and τ = 20.1182 fs (SnSe) were calculated from alignment with experimental data (see Methods for details). Already at low temperatures, the transport anisotropy along different directions is apparent, with the Y direction (p-type) matching the X direction (n-type) in SnS, while in SnSe, the latter is larger (Figure 6). The similarity of the plots follows from the similarity in the band structure, with a narrower band gap and  For SnS, the relaxation times at 500 and 750 K are set to 12.3156 fs and 4.7640 fs, respectively. For SnSe, the values at the same temperatures are τ = 7.6324 and 5.7991 fs, respectively. Progression from 500 to 750 K promotes the dominance of n-type over p-type, in agreement with other calculations. 25 The growth of the n-type maximum is noticeably more asymmetric in SnSe than that in SnS ( Figure  7), confirming the important role of interlayer spacing and structural anisotropy. In comparison, the electronic-only thermoelectric figure of merit is less anisotropic (see Figure

■ COMPARISON
Despite being structurally similar, Pnma-SnS is thermoelectrically inferior to Pnma-SnSe, due to higher thermal conductivity. As temperature increases, however, SnS becomes more similar to SnSe, with a steeper decrease in thermal lattice conductivity compared to SnSe (Figure 4). To better understand this effect, we have included a comparison with the high-temperature Cmcm-SnS phase. The latter results from compressing open zig-zag SnS layers ("accordion" geometry), achieving an additional in-layer Sn−S contact, accompanied by alternate layer antiparallel displacements. Its structural motif is not only relevant above the transition temperature, but it also plays a role of a "hidden" motif at lower temperatures through phonon modes that compress the accordion layers. We have therefore extended the comparison over the same temperature range as for the lower symmetry phases to deepen our understanding of how the structure−property relationship in this class of materials favors thermoelectricity. Thermal lattice calculations, geometry, and electronic structure for Cmcm-SnS are given in Section S5 and Figures S5−S8 in the SI. Lattice Thermal Conductivity. The direction-resolved lattice thermal conductivities are shown in Figure 8, separated into the cross-layer direction and the averaged intralayer direction. This highlights the difference between Pnma-SnS and Pnma-SnSe, with Pnma-SnSe less than half of Pnma-SnS across the temperature range. The Cmcm-SnS lattice thermal conductivity is closer to that of Pnma-SnSe than that of Pnma-SnS, indicating that structural similarities to Cmcm at higher temperatures are the reason for a better performance of Pnma-SnS. Anisotropy in κ latt is due to different contributions to κ latt along a, b, and c, particularly at lower frequencies ( Figure S9). In Cmcm-SnS ( Figure S9c) as in Pnma-SnS ( Figure S9a), the contributions to κ latt along X, and Y, Z, respectively, remain similar, while for SnSe, the contribution along Z is lower ( Figure S9b

■ DISCUSSION
The ranking of calculated thermal conductivities for Pnma-SnS and -SnSe is in good agreement with experimental results. Inclusion of three-phonon scattering and its effect on lifetimes ( Figure 9) seem therefore to account for the relevant physics, while isotope effects can be neglected as discussed elsewhere. 23 This is also in agreement with frequency lowering due to the larger Se mass and Sn−Se chemical bond softening. 25 Lower frequency phonon modes are major heat carriers as illustrated by the contribution to lattice thermal conductivity by frequency, with acoustic phonons becoming dominant from SnS to SnSe, as noticed above. Interestingly, the Cmcm-SnS phase appears closer to SnSe, pinpointing the central role of inlayer dynamics in reducing finite lifetimes, especially of lower optical frequencies.
Pnma-SnSe is characterized by lifetimes in the low-frequency region that are lower than Pnma-SnS, also reflected by shorter mean free path (MFP) contributions in Pnma-SnSe compared to Pnma-SnS. Compactification of MFP contributions to lower frequencies is also present in Cmcm-SnS, which, however, remains above SnSe values, due to larger group velocities in the Cmcm-SnS phase, in particular along the layering a direction. This entails the largest difference in lattice thermal conductivity between Pnma-SnSe and Cmcm-SnS. We argue that this is a consequence of symmetry and that temperatures just below phase transition would allow for group velocity rescaling, an effect that could be verified with lattice dynamics calculations. From the comparison of lattice thermal conductivities, phonon lifetimes, and mean free paths, phonon scattering near and across the transformation from Pnma to Cmcm appears key to thermoelectricity in this class of materials, while high symmetry seems to enforce higher phonon group velocities. Accordingly, Pnma-SnS turns into a better thermoelectric, as it becomes similar to Cmcm on increasing the temperature, as shown in this work. The increased importance of Cmcm at temperatures below the phase transition (T = 500 K) is illustrated in Figure S16 by isothermal−isobaric molecular dynamics simulations. In doing so, SnS becomes at the same time more like SnSe, thanks to in-layer deformation phonon softening.
Besides low lattice thermal conductivity, thermoelectric materials must be robust in their "electronic suitability." While Pnma-SnS and -SnSe are electronically closer than their lattice conductivities differ ( Figures S11 and S12), the former is superior in its Seebeck coefficient, while electrical conductivity is larger for SnSe. The a-axis power factors across the chemical potential range are very similar, whereas in-layer (b-and c-axis) values are larger for Pnma-SnS. Therefore, the main reason for the increase in the overall figure of merit of Pnma-SnSe is related to a larger decrease in lattice thermal conductivity. For Cmcm-SnS ( Figure S13), the low Seebeck is balanced by comparatively high electronic conductivities. Even if the interlayer (a-axis) power factor favors n-type and in-layer (band c-axis) favors p-type, the difference is less pronounced than that for the Pnma phase. This confirms that the role of Cmcm in lowering thermal transport lies mainly in its lattice dynamical influence, with electronics only moderately affecting its thermoelectric properties.
The maximal figure of merits ZTs for SnS and SnSe entail an increase in the layered direction as an n-type, whereas most experimental measures focus on in-layer p-type thermoelectrics, with values similar to the p-type figure of merits calculated in this paper. As the n-type maxima are obtained at larger chemical potentials, this would imply a further increase in the level of doping in real materials as a route to maximizing ZT beyond current reports.

■ CONCLUSIONS
In this work, we have calculated thermal end electronic transport for SnS and SnSe, using the frozen phonon and single relaxation time approximations. Phonons were calculated based on a transferable tight-binding (xTB 41 ) potential, which resulted in a reliable comparison with experiments, including the relative ranking of SnS and SnSe, anisotropy, and absolute lattice conductivity values. This improvement stems from larger supercells and from the xTB potential correctly capturing frequency softening at higher temperatures. Full figures of merit were obtained by including the electronic part via the Boltzmann transport equation based on the Wannier function interpolation of the electronic bands and single relaxation time, showing asymmetry between p-type and n-type regimes.
The lower thermal conductivity of SnSe compared to SnS is related to shorter lifetimes of optical modes, besides softer acoustic frequencies in SnSe. The inclusion of a complete analysis of the Cmcm high-temperature SnS phase suggests that this softening should be associated with in-layer b, c moduli compression, which indicates a dynamical role of a hidden Cmcm phase already at temperatures below phase transition. Once formed, Cmcm-SnS shows larger group velocities in the in-layer components, an intrinsic effect that contributes to preventing lattice thermal conductivity to reach even lower values.
Taken together, our findings provide a reliable, efficient, and scalable scheme to calculate thermal conductivity in tin chalcogenides, one that can be transferred to larger unit cells for composite material design and nanoinclusions. In combination with electronic transport calculations based on Wannier functions, a robust scheme for ZT thermoelectric rankings can be achieved, broadly suitable for thermoelectric material characterization.
Our analysis suggests that Pnma-SnS does become a better thermoelectric at higher temperatures, as it resembles more and more Cmcm-SnS. Dynamical in-layer compression is indicated as a key mechanistic element beyond thermal conductivity lowering, already before the Cmcm-SnS phase fully locks in. Leveraging this effect at lower temperatures is therefore key to engineering better thermoelectric materials based, here, on tin chalcogenides.
SnS and SnSe tight-binding (GFN-xTB) relaxed geometries; lattice thermal conductivity (pristine GFN-xTB forces); SnS phonons and comparison between QE and CP2K/xT; electronic figure of merit; full characterization of Cmcm-SnS (structure optimization and bands, phonons, thermal conductivity and phonon lifetimes, full figure of merit); comparison among SnS, SnSe, and Cmcm-SnS (contribution to κ latt as a function of frequency, ZT as a function of chemical potential, theoretical operational efficiency (TOE), and electronic transport calculations (Boltzwann)); comparison of RTA and LBTE results; convergence tests (lattice thermal conductivity and lattice thermal conductivities as a function of q-mesh size); and SnS: molecular dynamics simulations (PDF) programme for West Wales and the Valleys. This work used the Isambard 2 UK National Tier-2 HPC Service (http://gw4. ac.uk/isambard/) operated by GW4 and the UK Met Office and funded by EPSRC (EP/T022078/1). The authors acknowledge the support of the Supercomputing Wales project (ARCCA), which is part-funded by the European Regional Development Fund (ERDF) via the Welsh Government.